Bifurcation analysis of delay-induced resonances 
of the El-Nino Southern Oscillation 



By Bernd Krauskopf^ and Jan Sieber^, 

^Department of Engineering Mathematics, University of Bristol, UK, 
^Department of Mathematics, University of Portsmouth, UK 

Models of global climate phenomena of low to intermediate complexity are very useful 
for providing an understanding at a conceptual level. An important aspect of such models is 
the presence of a number of feedback loops that feature considerable delay times, usually 
due to the time it takes to transport energy (for example, in the form of hot/cold air or water) 
around the globe. 

In this paper we demonstrate how one can performed a bifurcation analysis of the 
behaviour of a periodically forced delay differential equation (DDE) in dependence on key 
parameters. As a concrete example we consider the El-Nino Southern Oscillation (ENSO), 
which is a sea surface temperature oscillation on a multi-year scale in the basin of the Pacific 
Ocean. One can think of ENSO as being generated by an interplay between two feedback 
effects, one positive and one negative, which act only after some delay that is determined by 
the speed of transport of sea-surface temperature anomalies across the Pacific. We perform 
here a case study of a simple delay-induced oscillator model for ENSO (introduced by 
Tziperman et al, J. CUmate 11 (1998)), which incorporates the two feedback effects and is 
parametrically forced by annual variation. More specifically, we use numerical bifurcation 
analysis tools to explore directly regions of delay-induced resonances and other stabihty 
boundaries in this DDE model for ENSO. 
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1. Introduction 

In climate science one finds a hierarchy of increasingly complex mathematical models, 
all the way from low-dimensional conceptual models in the form of ordinary differential 
equations (ODE) to global climate models with many millions of unknowns. Our interest 
here is in models of a low to intermediate complexity, in particular, those that feature 
delayed feedback mechanisms. The mathematical description of such a system, hence, takes 
the form of delay differential equations (DDEs). This type of model arises in other fields 
of science, for example, in mechanical engineering, in laser physics and in mathematical 
biology, where one finds communication delays (between optical elements or cells) as well 
as processing delays; see, for example, [1, 2]. In climate science, on the other hand, delays 
arise mainly in feedback loops due to the transport of mass or energy from one location on 
the planet to another. There are models that describe such transport phenomena directly in 
the form of partial differential equations [3]. However, in some situations it is advantageous 
to take a more conceptual point of view by considering only the resulting feedback effect 
itself, which is then subject to a delay due to the associated transport time. 

In this contribution we perform a case study of such a DDE model to demonstrate how 
tools from bifvu^cation theory can be brought to bear to tuiderstand its behaviotff. More 
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specifically, we study delay-induced resonances in the classical ENSO problem. This is an 
ideal showcase problem for bifurcation analysis because the ocean-atmosphere system of 
the Pacific has an oscillatory instability at the linear level due to reasonably-well identified 
feedback effects. The resulting oscillator interacts with the annual forcing to produce a rich 
set of resonances [4, 5, 6] in a wide region of parameters; see also [7]. Extensive research 
has created a hierarchy of models that isolate the processes affecting ENSO. We explore here 
one of the simplest of these models that still contains physically meaningful parameters, 
namely, the model introduced by Tziperman et al [8]. This model takes the form of a 
scalar DDE that lumps all coupled processes into two feedback terms, which act only after 
given delays due to the time it takes for thermocline anomalies (proportional to sea-siu^face 
temperature anomaUes) to travel across the Pacific. We provide a brief summary of the 
motivation behind this delayed oscillator modeld of ENSO in Section 2; for an account of 
the modelling of the processes involved at coarse to intermediate levels of complexity see 
the review [7] and the textbook [3]. 

Once a DDE model, such as that for ENSO, has been derived, the task is to understand 
its behaviour in dependence on relevant parameters. A quite straightforward method for 
investigating the eventual or attracting behaviour of a given mathematical model is its 
simulation by numerical integration. However, when one wants to find all possible behaviour 
by simulation, one is faced with a number of challenges. One generally needs to perform 
many simulation runs for different parameter values, for example, chosen from some 
sufficiently fine grid; each such run needs to be performed until transients die out sufficiently 
and care needs to be taken in the choice of the initial condition, especially in the presence of 
multistability. Note that for DDEs one needs an entire history segment as initial condition 
[9], making the search space rather large. Furthermore, there is the difficulaty of representing 
and identifying, ideally automatically, what the eventual dyanamics actually is. 

An alternative to simulation is bifurcation analysis. It was orginally developed as a 
mathematical tool to classify the long-term behaviour of nonlinear dynamical systems that 
are modelled by low-dimensional ordinary differential equations (ODEs) or maps. The 
underlying idea is that one studies equilibria and periodic orbits and their stability and 
other properties as a function of relevant system parameters. When a parameter is varied, 
one may observe well-defined qualitative changes of the dynamics, which are referred to 
as bifurcations. The goal of bifurcation analysis is to find the bifurcation diagram, which 
divides the parameter space into regions of qualitatively different behaviour. A bifurcation 
diagram is effectively a 'road map' of where in parameter space one may find which kind 
of behaviour, and what changes occur in the transition from one type of behaviour to 
another. In spite of their seeming simplicity, ODE models of quite low dimension (such as 
the well-known Lorenz and Rossler systems) can already give rise to intricate bifurcation 
diagrams with regions of stationary, periodic, quasiperiodic and even chaotic behaviour; 
see, for example, the textbooks [10, 11, 12] as entry points to bifurcation theory. 

Today, the bifurcation analysis of a given ODE can be performed routinely with freely 
available numerical tools, such as the packages AUTO [13] and MatCont [14]; see also the 
survey [15]. These tools only require knowledge of the right-hand side of the ODE and allow 
the user to perform a systematic exploration of the bifurcation diagram. This is achieved by 
finding and then tracking (or continuing) equilibria, periodic orbits, and their bifurcations to 
yield stability boundaries of the different types of solutions. The usefulness of bifurcation 
theory has been demonstrated in numerous studies of ODE models arising in diverse fields 
of science; see, for example, [16]. 

More recentiy, numercial bifurcation tools have been implemented also for more general 
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classes of systems. For complex models, which are typically based on partial differential 
equation formulations, numerical routines are available as libraries [17], and they have been 
used in studies of climate problems [3, 18]. Specifically for DDEs, numerical methods for 
their bifurcation analysis are today avalailable in the form of the packages DDE-Bif tool 
[19] and knut [20]. These numerical tools apply to general systems of DDEs, including 
DDEs with several delays and periodic forcing. In the same way as for ODEs, these tools can 
determine and track the stabiUty properties and bifurcations of equihbria and periodic orbits. 
As we will demonstrate, this means that numerical bifvffcation studies can be performed 
efficiently today also for models that feature delays; see also the survey [21]. 

The most obvious application areas of bifurcation analysis are those where systems 
under consideration have tunable parameters that the practitioner can adjust to achieve the 
desired behaviour [10]; concrete such examples are engineering and laser physics. The 
situation in climate science, on the other hand, is very different in that tuning a parameter 
is typically not possible in the real system. In spite of this fact, we argue that systematic 
parameter studies are still relevant in this context, namely to determine the sensitivity 
of the model under consideration to changes of parameters that are known only in an 
order-of-magnitude sense. The fact that the observed behaviour of a nonUnear system may 
depend critically on all parameters motivates the bifurcation analysis with respect to selected 
difficult-to-determine parameters; see [22] for an exemplary study. 

For our specific example, two parameters that affect the ENSO resonances strongly 
are the mean and the annual variation of the ocean-atmosphere coupling. The values of 
these two parameters are hard to ascertain in the real system; furthermore, they are also 
difficult to compare and convert between models of different complexity. Hence, it makes 
perfect sense to locate a chosen resonance tongue in the given model to restrict the values 
of these two parameters. When one attempts to locate a particular resonance one encounters 
the problem that the only places where one knows the frequency ratio analytically (or at 
least numerically without a scan of the complete parameter plane) is the small-forcing 
or the small-oscillation-amplitude regime. In these two regimes resonance tongues are 
extremely narrow, and convergence to the resonant orbits (if they are stable) is very slow. 
This makes resonances hard to observe and track with simulations in the small-amplitude 
regime (corresponding to the classical devil's staircase scenario [4]). In this paper we show 
how this difficulty can be overcome with numerical bifurcation analysis. More specifically, 
periodic orbits are tracked as solutions of a boundary- value problem (BVP) and, as such, 
these computations are unaffected by the instabiUty (or weak stabihty) of the periodic orbit 
under consideration or by its sensitive dependence on parameters. With this BVP approach 
we have tracked the relevant resonance tongues directly from their root point near the linear 
regime over a considerable region of the parameter plane. This bifmcation diagram, shown 
in Figure 6 of Section 4 below, constitutes the main result of our demonstration of the 
capabiUties of numerical bifurcation analysis in the context of cUmate science. 



2. A parametrically forced delayed oscillator model for ENSO 

We consider a conceptual model for the ENSO mechanism that simplifies the circulation of 
ocean surface water waves and the ocean-atmosphere interaction to a scalar DDE, which is 
periodically forced due to the annual variation of the strength of ocean-atmosphere coupling. 
We use the parameters and follow the notation of [8], where the dependent variable is the 
height anomaly h{t) of the thermocUne at the Eastern boundary of the Pacific. The height 
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5.75 (months) 


time it takes cold off-equatorial SST petur- 
bations from central Pacific to reach eastern 
boundary (first westward via Rossby wave, 
then eastward via reflected Kelvin wave). 


b 


1/120 (per day, 0.2535 per month) 


amplification factor of warm perturbations. 


c 


1/160 (per day, 0.1901 per month) 


attenuation factor of cold perturbations. 
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1/190 (per day, 0.1601 per month) 


attenuation due to dissipation. 


G{Kh) 


defined in Eq. (2.2) 


ocean-atmosphere coupling (saturation non- 
linearity). 


K 


K= Ko+AKSin{cot) 


annually varying slope of G at /! = 0, 


(0 


2k /\2 (per month) 


frequency of annual variation (period T = 12 
months). 


b± 


b+ = -i,b-=-\ 


maximum and minimum of saturation nonlin- 

earity in G. 




Table 1: Values and interpretation of parameters in (2.1)-(2.3) 


anomaly h (in dimensionless units) satisfies 






= bG[K{t-Xi)h{t-Xi)\-cG[K{t-X2)h{t-T2)]-dh{t), where (2.1) 


G{x) 


{b+tw\\{x/b+) ifx>0, 
1 &_tanh(x/fo^) ifx<0, 


(2.2) 


Kit) 


= kQ + dkSin{(Ot). 


(2.3) 



The right-hand side of (2.1) combines positive feedback due to the warm equatorial sea- 
surface temperature (SST) perturbation in the central Pacific, negative feedback due to the 
cold off-equatorial SST perturbation in the central Pacific, and dissipation. The parameters 
and their values used here are from Tziperman et al [8], and they are listed in Table 1. 

The idea that a dynamical instability is the basic mechanism behind ENSO was first 
proposed by Bjerknes in 1969 [23]. The review by Neehn et al [7] describes the hierarchy of 
models that have been developed during the following decades. The more complex models 
in this hierarchy help to establish which physical processes are connected. The simpler 
models, such as (2.1) considered here, take these connections for granted and treat them 
as lumped positive or negative feedback terms in the right-hand side. Figure 1 sketches 
the physical processes that enter the right-hand side of (2.1) and their connections. The 
height of the thermocline is a proxy for the sea-surface temperatvffe. In equilibrium the 
thermocline is low at the eastem boundary of the Pacific and high at the western end. This 
difference maintains an atmospheric convective loop above the equatorial ocean, which 
points westward at the ocean surface because the relatively cold water at the eastem end 
cools the atmosphere, causing air to sink (and vice-versa, the relatively warm water at the 
western heats up the atmosphere, causing air to rise). The dependent variable h in (2.1) 
measures the anomaly of the thermocline at the eastem boundary, that is, the deviation from 
equilibrium at this point. 
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Figure 1: Schematic of tiie two feedbaclc meclianisms that motivate and are modelled by the ENSO model (2.1). 
The .shading of the ocean corresponds to the thermocline (and sea-surface temperature) in equilibrium: darker 
means lower (colder). Equation (2.1) only describes the anomaly of the thermocline, h, at the eastern boundary. 
The labels /jec and /Jqc refer to the thermocline anomaly in the equatorial and the off-equatorial central pacific, 
respectively, which feed into the feedback mechanisms. Regular aiTows indicate same-sign re-inforcement, bold 
arrows oppposite-sign re-inforcement, and the label on each arrow show the time delay of the interaction. 




Figure 2: Profile of coupling function G(Kh) for IC G 0, . . . , 5.5, which are the values that are encountered in the 
parameter range {dt.ko) G [0,2.5] x [0,3]. 

A disturbance in the atmospheric convective loop, that is, a slowing down of the 
westward surface wind (an eastward wind stress of the ocean surface) in the central Pacific 
causes surface water transport toward the equator Surface water is warm, causing an increase 
of h in the equatorial central Pacific, and a decrease of h in the off-equatorial central Pacific. 
The positive thermocline disturbance at the equator travels eastward along the equator (in a 
Kelvin wave), and hits the eastern boundary after time Ti . The corresponding increase in 
sea-surface temperature at the eastern boundary slows down the atmospheric convection 
(the relatively cold water at the eastern boundary is a driver of the convection). This creates 
the first, positive feedback term in (2.1) with the delay Ti (note that b > 0). 

The off-equatorial depression of h travels westward and toward the equator (in a Rossby 
wave), reflects at the western Pacific boundary, and then travels as a Kelvin wave eastward, 
reaching the eastern boundary after time T2. The depression of h corresponds to a cooling 
down of the sea surface such that atmospheric convection is enhanced, creating the negative 
delayed feedback term in (2.1) with the delay T2 (note that c > 0). 

The above effects work all at the linear level of the disturbance of the equilibrium. 
The nonlinearity in (2.1) comes from saturation and asymmetry in the ocean-atmosphere 
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coupling instability, as expressed by the function G in (2.2) that is shown in Figure 2. 
Moreover, the slope of G at = varies throughout the year, being strongest in the northern 
summer and weakest in the winter; see [8] for a realistic time profile of the annual variation. 
This annual variation is expressed in the sinusoidal variation (2.3) of the pre-factor k in the 
argument of G. Since only the ampUfication factor varies in time, the oscillator described by 
(2.1) is parametrically excited; in particular, the equilibrium h = Qof the unforced system 
(for dk = 0) corresponds to a trivial periodic solution in the periodically forced system (for 
dk > 0). 

3. Delay differential equations and their bifurcations analysis 

This section is a short review of the general theory of dynamical systems with delays as 
needed for the bifurcation analysis of the ENSO model in Section 4, which the reader 
may wish to consult already for illustrations of the concepts introduced now. To keep the 
exposition specific to the ENSO model under consideration, we present the theory for the 
case that the right-hand side of the DDE depends on the current value of the dependent 
variable, as well as on two values from the past, with delays Ti < T2. Furthermore, we 
allow the right-hand-side of the DDE to depend explicitely on time, which accounts for the 
periodic forcing of (2.1). Hence, we consider a DDE of the form 

^x(f) =/(?,x(0,x(?-Tl),x(?-T2),Tj), (3.1) 

where x e K" is an n-dimensional state and r] e K™ is a multi-parameter; one also refers to 
K" as the physical space. The fimction 

/ : M X M" X R" X X M"* -!> R" 

describes the right-hand-side, and it is assumed to be sufficiently smooth. Note that (3.1) 
is non- autonomous because / depends explicitely on time t. We consider here only the 
case of periodic forcing, that is, /(f, 77) = f{t + Tf , t]) where Tf is the forcing 
period. An important special case is that the forcing ampUtude is zero, so that system (3.1) 
is autonomous. We assiune here that the forcing amplitude is one of the components of 
the parameter vector Tj, so that the limiting autonomous case can be foimd as a subset of 
parameter space. 

We proceed to discuss the basic properties of the DDE (3.1), both for the autonomous 
and the non- autonomous cases, which form the basis of the bifurcation analysis of the 
ENSO model in Section 4. We remark that the theory present here is also valid for more 
general DDEs, including those with distributed delays, as long as all delays are constant 
and bounded. Full technical proofs can be found in the classical textbooks by Hale and 
Verduyn-Lunel [9] and Diekmann et al [24]; the textbook by Stepan [1] is also rigorous, 
but provides a shorter route to those parts of the theory that are necessary for practical 
appUcations. 

Because of the presence of delayed terms, the DDE (3.1) does not define a dynamical 
system on the physical space M". Rather, its phase space is the infinite-dimensional space 
M X C([— T2,0];M"); here E represents time and C([— T2,0];M") is the space of continuous 
functions over the maximal delay interval [— TmaxjO] = [— T2,0] with values in the physical 
space W. This means that in order for the initial- value problem to be well defined one needs 
to prescribe an entire fimction : [— T2,0] — )• K" as the initial condition. Similarly, after 
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solving (3.1) from time to time / the current state is the function qt : [— T2,0] — >^ M" given 
by qt{s) =x{t-\-s), s [— T2,0]. These trajectories qt in the infinite-dimensional phase space 
M X C([— T2,0];IR") depend smoothly on the initial condition q^,, meaning that classical 
dynamical systems theory can be appHed to the DDE (3.1) in M x C([— T2,0];M"). When the 
forcing amplitude is zero the time axis decouples and it is sufficient to study the autonomous 
DDE with phase space C([-T2,0];M"). 

(a) Equilibria and their bifurcations 

We first consider the case that / is autonomous, which is of interest on its own as well 
as in the context of periodic forcing where it corresponds to zero forcing amplitude. An 
equilibrium or steady state is an initial fvmction go : [—^2,0] R" in phase space such 
that qt{s) = qo{s) for all s G [—^2,0] and t >0. Hence, qo{s) equals a constant xq for all 
■s € [— T2, 0] and the constant vector xq G W can be determined from the equation 

f{xo,xo,xo,-n)=0, (3.2) 

which, as for ODEs, is a system of n algebraic equations. This means that equilibria of 
DDEs can be foimd and tracked in a parameter as roots of (3.2) by standard continuation 
methods. 

To determine the stability of an equilibrium qo one moves it to the origin and linearises 
the DDE. The result is a Unear DDE of the form 

x{t) = f{t,x{t),x{t-Xi),x{t-X2),'r]) =Ax{t)+Bx{t-Ti)-\-Cx{t-X2), 

where A, B and C are n x n matrices. The stability of the origin is determined by the spectrum 
of the eigenvalues of the linearisation, which can be found as the roots of the characteristic 
function 

;t;(A) = det(A/-A-Be-^^i +Ce-^'^). (3.3) 

The characteristic function ;|f of a linear autonomous DDE (of the type considerd here) has 
at most finitely many roots with non-negative real part [9, 1]. This means that the stability 
theory of equilibria of DDEs is very similar to that for ODEs. An equilibrium is stable if 
all eigenvalues of its linearisation have negative real part. It changes its stability type when 
eigenvalues cross the imaginary axis of the complex plane. When a single parameter is 
varied there are two typical scenarios, which give rise to the standard local bifurcations of 
equiUbria of codimension-one [10, 11, 12]: 

• when a single real eigenvalue goes through zero the bifurcation then one encounters, 

in the generic case, a saddle-node bifurcation, where two equilibria meet and disapper 
(or are created); in the presence of additional properties of the equation (such as 
symmetries) one encounters a branch bifurcation, where two additional equilibria 
bifurcate; 

• when a complex conjugate pair of eigenvalues moves across the imaginary axis of the 
complex plane then one encounters a Hopf bifurcation from which a small periodic 
orbit bifurcates. 

These bifurcations of autonomous DDEs can be detected and tracked as for ODEs by 
fixing the real parts of roots of the characteristic function x in (3.3) to zero and freeing 
a second parameter (such that tj e M?). In this way, one can compute curves of saddle- 
node/pitchfork and Hopf bifurcations in a two-dimensional parameter space. The curves 
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divide the parameter plane into regions where different numbers and types of equilibria 
exits. 

{b) Periodic orbits and their bifurcations 

An initial function qo is a periodic point of the DDE (3.1) if it satisfies qTri^) = for all 
•s G [— T2,0] and some Ty > 0, where Tr is called the (minimal) period of qo. Practically, one 
finds a periodic point of (3.1) as a periodic orbit, which is a closed curve F : [0, 1] K" in 
the physical space satisfying the periodic boundary-value problem 

t{t) = 7r/(/,r(0,r(?-Ti/7r)„od[0,il,r(r-T2/7r)^od[0,i].n) for; e (0,1], 

r(o) = r(i). 

The notation r(f — T,7rr)niod[o,i] means that for combinations of /, Tj and Tr for which 
t — Tj/Tr is outside of the interval [0, 1] we choose the appropriate integer v such that 
t - Tj/Tr + V e [0, 1] and evaluate r(f - Tj/Tr + v). From any solution F of (3.4) one can 
find periodic points q of period Tp through the relationship 

q{s) = r(i/7r)mod[o,i] ^ ^ [—^^2,0]. (3.5) 

When the DDE is autonomous then the period Tp in (3.4) is one of the unknowns that must 
be found together with F. This can be done in the same way as for ODEs by the adding 
one scalar equation, a so-called phase condition [25], which determines Jp and makes the 
solution r of (3.4) locally unique. Good starting points for periodic orbits of autonomous 
DDEs are Hopf bifurcations where small-ampitude periodic orbits branch off. Their period 
Tr is determined by the imaginary part of the pair of purely complex conjugate eigenvalues 
in the Hopf bifurcation. 

The situation is different when the DDE is periodically forced. Then there are no 
equilibria, and the simplest invariant object is a periodic orbit F whose period equals the 
forcing period Tf, that is, Tr = Tf is known in (3.4). In fact, each equilibrium qo{s) = 
XQ, s G [— T2,0] of the autonomous system for zero forcing amplitude gives rise to a trivial 
periodic orbit of the form Fo{t)=xo. The trivial periodic orbit Fq can then be continued 
as a periodic orbit F of the periodically-forced DDE when the forcing is increased from 
0. For the special case that the forcing is parametric, as in the ENSO model (2.1), one 
finds that F{t) = Fo{t) = remains constant for all forcing amplitudes; however, in general, 
r(;) is a non-constant periodic function when the forcing is nonzero; note that generically 
the period of any periodic orbit of a periodically-forced DDE is an integer multiple of the 
forcing period. Periodic orbits can generally not be found analytically, but (3.4) can be 
solved nvunerically with the same numerical methods as one computes periodic orbits of 
ODEs [21]. 

The stabihty of a periodic orbit F is determined by its Floquet multiphers Hi, which are 
the eigenvalues of the time-Tp map of the point q on the periodic orbit F given by (3.5). If all 
jii are inside the unit circle of the complex plane then F is stable, and the eigendirections of 
the Floquet multiphers outside the unit circle correspond to unstable directions. Similarly to 
the case of equiUbria, for a time-periodic DDEs one can construct a characteristic function 
XpilJ-) such that the Floquet multipliers of F are the roots of Xp [20, 26]. Generally, the 
characteristic function Xp and its roots are found numerically, after one has computed F 
and its period Tp. For a DDE there are infinitely many roots of Xp and, hence, Floquet 
multiphers Hi, but they are also discrete and have the origin of the complex plane as their 
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only accumulation point; in particular, there are at most finitely many unstable Floquet 
multipliers. As for ODEs, a periodic orbit changes its stability when Floquet multipliers 
crosss the unit circle, and this gives rise to the same well-known local bifurcations of 
periodic orbits. When a single parameter is varied there are three typical scenarios, which 
give rise to the standard generic local bifurcations of periodic orbits of codimension-one 
[10, 11, 12]: 

• when a single real Floquet multiplier goes through +1 then one encounters a saddle- 
node bifurcation of limit cycles, where two periodic orbits meet and disapper (or are 
created); 

• when a single real Floquet multipUer goes through — 1 then one encounters a period- 
doubling bifurcation, where a periodic orbit of twice the period bifurcates; 

• when a complex conjugate pair of eigenvalues goes through the unit circle at e^'^^"' 
then one encounters a Neimark-Sacker or torus bifurcation, where a smooth in- 
variant torus bifurcates (provided that low-order resonances are avoided, that is, 
a 1 i 1 2 i 3x 

' '2'3'3'4'4'^' 

These bifurcations of DDEs can be detected and tracked as for ODEs by appending the 
above conditions on the Floquet multipliers to the continuation of the periodic orbit. In 
this way, one can compute curves of saddle-node of limit cycle, period-doubUng and torus 
bifurcations in a two-dimensional parameter space. 

(c) Resonances and quasiperiodicity 

Resonant periodic orbits are special periodic orbits that correspond to the locking of two 
different oscillations of the system. They are organised in DDEs in resonance tongues and 
resonance surfaces in the same way as in ODEs; see already Figure 6 and, for example, 
[27, 28 1 for the general theory. Along a torus bifurcation curve in a parameter plane a 
Floquet multipher g^^aai j^oves along the unit circle, that is, the quantity a changes. When 
a = A:/^ e Q a two-parameter family of A: : £ resonant periodic orbits (a so-called resonance 
surface) branches off. The projection of this surface into the parameter plane is a resonance 
tongue, which is a region in the parameter plane where periodic orbits exist that are A: : £ 
locked to the original periodic orbit (or the forcing period in a forced DDE such as the 
parametrically forced ENSO model). Inside the resonance tongue there are two periodic 
orbits, one of which is always unstable. The root point of the resonance tongue is attached 
to the bifurcation curve. Close to the root point the two periodic orbits form k : i torus knots. 
In a cross section transverse to the flow one finds a period-^ orbit; the return map to the 
cross section acts on it by mapping each point to its kCa neighbour in the counter-clockwise 
direction. If the denominator i is also greater than 4 both periodic orbits lie on an invariant 
torus (close to the root point this is guaranteed by general theory [12]), one periodic orbit is 
stable within the torus one unstable. 

Similarly, one finds resonance tongues that are rooted at a line in the parameter plane 
where the forcing is zero. At such a root point of the k : I resonance tongue the unforced 
system has a periodic orbit F such that the ratio of its period and the forcing frequency 
7y is a = kji e Q. 

By contrast, when a is irrational, that is, a ^ Q, then there is a smooth curve in the 
parameter plane starting at the corresponding point on the torus bifurcation curve, or line of 
zero forcing, along which one finds a torus that is densely filled by a quasiperiodic trajectory. 
The overall scenario is one of a complicated interplay between locked and quasiperiodic 
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dynamics on the bifurcating torus. In particular, if one plots the rotation number along 
any parameter path that intersects the resonance tongues one finds the well-known devil's 
staircase consisting of infinitely many intervals where the rotation number is constant (inside 
the resonance tongues) as also observed in simulations of ENSO models [4, 5, 6]. 

A resonance tongue is bounded near its root point by saddle-node of limit cycle bifurca- 
tion curves of the two orbits existing inside. An alternative method to the simulations as 
done in [4, 5, 6] is to compute the : f resonance surface by using the periodic BVP (3.4) 
with period Tp = tTf and two free parameters. In this way one obtains a surface consisting 
of all locked periodic orbits inside the resonance tongue. This compuation can be started 
at any point with rational a on the torus bifurcation curve or the line where the system 
is unforced. The continuation of the locked periodic orbits is numerically robust, and the 
resonance tongue is obtained by projection onto the parameter plane. This method has been 
introduced in [28, 29] for the computation of resonance surfaces in maps and forced ODEs; 
Appendix A explains how it can be adapted and implemented for a general DDE of the form 
(3.1). 



{d) Ready-to-use software for the bifurcation analysis of DDEs 

Finding and continuing equihbria and periodic orbits and their codimension-one bifurcations 
has been implemented in the saftware packages DDE-Bif tool [19] and knut [20]; see 
also the survey [21]. The basic underlying idea is to employ very similar techniques as 
for the well-established numerical bifurcation analysis of ODEs. To this end, the infinite 
dimensional DDE needs to be discretised to find a relevant number of right-most eigenvalues 
of equihbria and largest Floquet multiphers of periodic orbits, respectively. The boundary 
value problem (3.4) defining periodic orbits is solved by both tools with Gauss collocation 
as in the numerical bifurcation analysis package AUTO [25]. 

DDE-Bif tool is a collection of Matlab functions that perform each necessary task 
and that can be called by the user. For example, there exist functions for calculating the 
eigenvalues of a given equilibrium, the Floquet multipliers of a periodic orbit, for converting 
a periodic BVP into a large (discretised) nonhnear system, for solving nonhnear equations 
with Newton iteration, for adapting the mesh on a given periodic solution, and for tracking 
solutions of nonlinear systems if one has one more free parameter than equations. Curve 
tracking of equilibria in a single parameter, of codimension-one bifurcations of equihbria in 
two parameters, and of periodic orbits of autonomous DDEs in a single parameter have been 
implemented. Codimension-one bifurcations of periodic orbits can be detected, but tracking 
them requires one to define a suitable extended system manually. Because DDE-Bif tool 
permits one to add an arbitrary number of conditions and parameters to the system, it is 
possible to track solution surfaces instead of curves. We make use of this feature for the 
computation of resonance tongues; see Appendix A. 

The package knut (formerly known as PDDE-Cont) has implemented curve tracking 
of equilibria and their codimension-one bifurcations, and monitoring of their stability. For 
periodic orbits, it has implemented curve tracking (in particular, it permits one to branch 
off from Hopf bifurcations or from zero forcing), eigenvalue monitoring and detection and 
curve tracking of codimension-one bifurcations. The numerical methods employed in knut 
are summarised in the review article by Roose and Szalai [21] and in [20]. The user needs 
to specify only the right-hand side /, and can then control the curve tracking and branching 
by providing an initial guess and by specifying which parameters to vary. 

Article submitted to Phil. Trans, of the Royal Society 



Bifurcation analysis of delay-induced resonances ofENSO 



11 



Re A 




0.2 
0.1 

ImA 

-0.1 
-0.2 





^^increase kp 




li 

H 



-0.5 -0.4 -0.3 -0.2 -0.1 



0.1 0.2 

Re A 



Figure 3: Stability of the trivial equilibrium /! = and of bifurcating secondary equilibria of (2.1) for i4 = as a 
function of k^. Panel (a) shows the real parts of the eigenvalues of the equilibria as a function of k(j, and panel (b) 
shows the locus of the two leading eigenvalues in the complex plane; the eigenvalues of /i = are plotted in black 
and those of the bifurcating equilibria in grey. 



4. Bifurcation anlyisis of the ENSO model 

The main conclusion of Tziperman et al [8] was that in model (2.1) the peak of the ENSO 
is locked to the end of the year, which corresponds to the minimum of the excitation K(t). 
That is, if Kit) ^ko + dkSm{t7z/6) such that the period of (fis 12 (time is in unit of months), 
then the maximum of the thermocline anomaly, h, occurs at f sa 8 . . .9 + 12£ where i is the 
denominator of the resonance. In [8] the parameters were chosen such that £ = 3, so the 
ENSO peak occured every three years. The authors also stated that in all their simulations 
they found the locking to the peak to the end of the year to be extremely robust with respect 
to all parameters. 

We demonstrate here how the phenomenon of delay-induced resonances can be studied 
systematically by means of a bifurcation analysis of (2.1). More specifically, we determine 
allk : £ resonance tongues with ^ < 1 1 in the relevant region of the (ii^-,fco) -plane of mean 
Icq versus annual variation dk of the amplification factor in (2.3); they are chosen as the 
bifurcation parameters since the interaction between ocean and atmosphere is known only 
qualitatively as expressed by the function G shown in Figure 2. In contrast, confidence in 
the other parameters of (2. 1) is higher: the delays Ti and T2 are known with good accuracy, 
and the relation between damping and feedback effects can be estimated directly from 
observations (involving only oceanic quantities). These other parameters are fixed at the 
values from [8]; see Table 1. 



(a) The autonomous regime for dk = 

The starting point of our bifurcation analysis is the autonomous case where d^ — 0. Then 
/; = is an equilibrium irrespective of the value of Icq, and we know that it is stable for 
k() ~ 0. We proceed to continue this equilibrium with the package DDE-Bif tool [19] over 
the range ^0 G [0,3] while monitoring its stability. The behaviour of the eigenvalues of 
/z = is illustrated by the black curves in Figure 3. Panel (a) shows how the real part of 
the eigenvalues changes with Icq, and panel (b) how the two leading eigenvalues move 
in the complex plane with increasing ^o- The trivial equilibrium /i = is initially stable 
and then loses its stability at feo ~ 1 -426 in a Hopf bifurcation H, where a complex pair of 
eigenvalues crosses the imaginary axis. At ^0 ~ 2.158 the complex pair meets on the real line 
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Figure 4: Single-parameter bifurcation diagram (a) of the autonomous system (2.1) for df. = and varying Uq. The 
equilibrium ft = loses stability at the Hopf bifurcation H; a branch of secondary equilibria emerges from the 
pitchfork bifurcation PF. Panel (b) shows the T along the branch of stable periodic orbits that bifurcates from H, 
and panels (c)-(g) show time profiles of the indicated low-order resonant periodic orbits along the branch. 



of the complex plane and splits up into two real eigenvalues. One of these real eigenvalues 
decreases and crosses the origin of the complex plane at k(, m 2.526 (the point is labelled 
B in Figures 3 and 4). Physically, at B the positive feedback terms in (2.1) exactly cancel 
the damping and the negative feedback for /i = 0, and a branch of two nonzero equilibria 
bifurcates and exists for past B. We continued these equilibria, and their eigenvalues are 
shown as grey curves in Figure 3. There is a pair of positive real eigenvalues and, hence, the 
secondary equilibria are unstable in the range G [2.526,3]. 

Figure 3(a) shows the one-parameter bifurcation diagram for k{) e [0, 3]. The equilibrium 
h = loses its stability at H, and at the branch point B the branch of secondary, unstable 
equilibria emerges. Also shown in Figure 3(a) is the branch of bifurcating periodic orbits 
as computed by continuation with DDE-Bif tool from the Hopf bifurcation point H. The 
periodic orbits are represented by plotting their minina and maxima; they are all stable, as 
was checked with the computation of the Floquet multipliers (not shown). Of special interest 
are those periodic orbits along the branch that give rise to low-order resonances with the 
forcing period of Tf — 12 months (once the annual variation of the amplification is turned 
on, that is, dk is increased from 0). Figure 4(b) shows the period 7r of the periodic orbits 
along the branch, where all k: i resonant periodic orbits with i< II are indicated — from 
the 1 : 4 resonance very close to the Hopf bifurcation H to the 1 : 6 resonance at ko « 2.982. 
Panels (c)-(g) show the respective time profiles of these resonant periodic orbits. 

Figure 4 gives a complete picture of the autonomous dynamics of (2.1) for dk — 0. Notice 
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further from Figure 3(a) that all but the two leading eigenvalues have quite large negative 
real parts. This suggests that the dynamics of the DDE (2.1) behaves like a two-dimensional 
ODE with parametric periodic forcing. 



(b) The dynamics with periodic forcing for di^y Q 

We now perform a bifurcation analysis of the periodically-forced DDE (2.1) for dit_ > 0. The 
trivial periodic solution h = Q for dk — gives rise for c/jt > to a periodic orbit with period 
Tf. Owing to the parameteric nature of the periodic forcing, this periodic orbit has zero 
amplitude, that is, it remains the trivial solution h = 0. As was the case for the equilibria in 
the autonomous case, the main question is when /i = is stable. This can be determined by 
its continuation in d^ for fixed Icq. Figure 5 shows how the Floquet multipliers of the trivial 
periodic orbit h = change as d^ is increased for fixed fco = 1 .8, which is the value that was 
used in [8]. For [kQ.dk) — (1.8,0) the periodic orbit /i = is initially unstable, since a single 
pair of complex conjugate Floquet multipliers lies outside the unit circle. At d^ ~ 1.659 
this complex conjugate pair crosses the unit circle and /i = undergoes a torus bifurcation, 
labelled T, and becomes stable. When dk is increased further, the complex conjugate pair 
comes together on the real line of the complex plane sAdk~ 1.810 and then splits up into 
two real Floquet multipliers. Very soon thereafter, at d^ « 1.814, one of them decreases and 
crosses the unit circle at — 1 and, hence, h = Q undergoes a period-doubling bifurcation, 
labelled P, and becomes unstable again. 

We now consider the bifurcation diagram of (2. 1) in the (t/^., A:o)-plane, which is shown in 
Figure 6. The torus bifurcation and the period-doubling bifurcation have been continued as 
curves. The torus bifurcation curve T starts at the Hopf bifurcation point on the automonous- 
limit line dk = and it ends on the period-doubling bifurcation curve P in a 1 : 2 resonance. 
Below the curve T and to the left of the curve P the trivial periodic orbit /z = is stable. To 
the right of the curve P the periodic orbit /; = is unstable. 



(c) Resonance tongues in the ENSO model 

Above the torus bifurcation curve T one finds invariant tori and associated k : £ resonance 
tongues, and Figure 6 shows all resonance tongues for £ < 1 1. As we have already seen in 
Figure 4, along the line dk = and above the Hopf bifurcation point H we find five root 
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Figure 6: Two-parameter bifurcation diagram of tlie ENSO model (2.1) in the (dt,A:o)-plane {d^,kQ). All k : I 
resonance tongues with £ < 1 1 are shown. They are rooted at points on the autonomous boundary = and on 
the torus bifurcation curve T, which ends on a period-doubling curve at a 1 : 2 resonance. The trivial solution 
/! = is stable in the region below the curve T . The inset shows an enlargement near the tip of the 1 : 3 resonance 
tongue; the point coiresponding to the parameter values used by Tziperman et al [8] is marked by a cross. 

points of these resonance tongues. Additional root points of resonance tongues are found 
along the torus bifurcation curve T . Specifically, we observe that the rotation number a 
gradually increases until the torus bifurcation curve terminates in the 1 : 2 resonance on 
the period-doubling bifurcation curve P\ along T we find nine more root points of k : £ 
resonance tongues with ^ < 1 1 . The resonance tongues in Figure 6 have been computed as 
resonance surfaces with the method described in Appendix A; they were then projected onto 
the {dk,ko) -plane to yield the open regions of locked tori and their boundaries. Notice how, 
as theory states [27, 28], the resonance tongues become very narrow with the denominator 
£. We found that the corresponding k : £ locked periodic orbit of the ENSO model is stable 
and, hence, observable practically everywhere throughout the considered parameter region 
(exceptions are the parts of the narrow tongues 3:11,3:10, 4:9, 5: 11 in the area where 
dk > 2). The cross in Figure 6 shows the parameter point chosen by Tziperman et al [8]; it 
is slightly to the left of the 1 : 3 tongue (see the inset), which explains that these authors 
observed non-periodic motion closely resembling a period-three orbit. 

The other main feature observed in [8] is that the peak of the oscillations appeared to be 
locked to the end of the year, that is, the maximum of h{t) occurs at time f w 8 . . .9 + 12^. 
To investigate how general this observation is Figure 7 displays the profile of all stable 
periodic orbits in the main 1 : 3 and 1 : 4 resonance tongues, together with the forcing. We 
find clear evidence of phase locking to the forcing in both cases, where the stable periodic 
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Figure 7: Locking of ENSO peak to the end of the year for all stable periodic orbits in the main tongues of 1 : 3 (a) 
and 1 : 4 resonance (b), where grayscale represents the value of df. for which the periodic orbit is found. The black 
sinusoidal curve at the bottom indicates the phase of the forcing. 

orbit reach its maximum shortly before the forcing reaches its minimum. Since the delays Ti 
and T2 enter directly into the phase of the forcing of the ENSO model (2.1), it is clear that 
the two delays strongly influence the location of the maximum of h relative to the forcing 
phase. Figure 7 is evidence that the phase between h and the forcing is nearly independent 
of the parameters dk and ko. This observation provides a certain justification for a reduction 
of the ENSO model even further to contain only the delays and not the periodic forcing; see 
[30] for a demonstration of how resonances can be explored with Boolean delay equations. 



5. Conclusions 

We presented a bifurcation study of how delay-induced resonances organise the parameter 
plane of mean and variation of the annual forcing of a model of ENSO with delayed feedback 
and forcing. The main message is that, colloquially speaking, the bifurcation theory for 
the considered class of DDEs is exactly as that for ODEs, but embedded into an infinite 
dimensional phase space. Hence, equilibria, periodic orbits and their bifurcations can be 
detected and tracked in DDEs much as for ODEs, and software packages exist for this task. 
As our case study demonstrates, numerical bifurcation anlysis is today perfectly feasible for 
periodically forced DDEs, and it allows one to obtain insights beyond what is possible with 
long-time simulations. 

Bifurcation analysis of equilibria and periodic orbits can give evidence for the presence 
of deterministic chaos via the detection of bifurcations associated with routes to chaos. A 
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specific example is chaotic dynamics due to the break-up of tori when resonance tongues 
overlap. In the ENSO model the resonance tongues do not overlap in the parameter region 
we explored and, hence, no chaotic dynamics was found. However, when the forcing is large 
enough, resonance tongues will generally overlap, and the method we used of computating 
resonance tongues as surfaces is especially suitable in this case. Hence, chaotic dynamcis 
can be studied in DDE models in considerable detail. An issue in this regard is that — even 
for systems much simpler than climate models — it is often very hard, if not impossible, 
in practice to distinguish between deterministic chaos and noise-induced fluctuations in 
observed data. One reason for this may lie in a lack of hyperbolicity of the deterministic 
chaotic attractor, which makes it non-robust (that is, sensitive to arbitrarily small changes in 
the parameters); this non-robustness gets 'smeared out' if one adds random disturbances to 
the deterministic model; see [31] for a study relevant to climate dynamics. 
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Appendix A. Computation of resonance tongues as surfaces 

We now present technical details of how we compute a resonance surface of a DDE of the 
form (3.1); our method is based on the principle described in [28, 29] for ODEs and has 
been implemented on top of DDE-Bif tool. 

(i) Branching off from a resonant torus bifurcation point 

We consider ak: £ resonant point on a previously computed torus curve in a parameter 
plane, where the periodic orbit F undergoing the torus bifurcation has a Floquet multiplier 
^±2nai a = k/£ where <k < £ are coprime integers; we represent this point as the 
tuple {x,y,a,ri). Here t] is assiuned to be of dimension two, and the periodic orbit F is 
the solution of the periodic BVP (3.4). Let z{t) be the eigenfunction corresponding to the 
Floquet multiplier g^jna. ^^^^ ^ period £ (since k and £ do not have a common divisor). 
We now choose a small radius p, and consider as initial guesses for the initial circle (of 
radius p) of locked periodic orbits the family 

xo^^{t)=r{£t)+pcos(l)Rez{it)+psin(l)Iraz{£t), (j) e [0,2k/£], te [0,1], 

which is parametrised by the angle (j). Note that here we have rescaled time again so that 

xo^^j, has period 1. We now perform a Newton iteration where we keep the two-dimensional 
system parameter t] (initialised to its value at the resonance point) free but fix the integral 
scalar products 



Together with requiring that x,j, satisfies the periodic boundary- value problem (similar to 
(3.4) but with Tr = ^7/^) 



system (A 1)-(A 2) is an equation in for the variables (x^ , TJ ) with a locally unique solution 
for every <j) and (sufficiently small) initial radius p. 
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(ii) Branching off from an autonomous periodic orbit 

Creating the initial circle for resonance surfaces starting from an autonomous oscillation 
for zero forcing is slightly different from the case of branching off from a torus curve. 
Suppose that at 7] = (771 , 772) — (771 ,0) the period Tp of a periodic orbit Y and the forcing 
period Tf have a ratio k/£, where < A: < £ are again coprime; here we assume that the 
second component TJ2 of the parameter 77 is the amplitude of the forcing. Then we choose 
as our initial guesses xq.,;, (f ) = x(l{t + ((>/ {^k))) (again, rescaling time such that the period 
of xo,(^ equals 1). Two additional conditions are given by fixing the initial forcing 772 to the 
small radius p, and by fixing the phase of the initial solution x^: 

= / xo^^{tfx^{t)dt, 772= p. (A3) 

J 

Solving (A 2) combined with (A3) results in an initial topological circle of locked orbits 
(rji,rj2)), parametrised by (j) e [0,2n/r]. 

(iii) Growing the surface circle by circle 

Once we have an initial circle of solutions x^, on the resonance surface we can continue 
the surface by computing nearby circles (which are generally only topological circles). We 
parametrise the surface locally by the angle on the circle and by its orthogonal complement 
in the tangent space of the solution surface. An equidistributed mesh of solutions on an 
already computed topological circle is given. We pick a point {x,r\)ou at angle ^ on this 
circle. The tangent space in (x, 77)oid to the resonance surface is two-dimensional, and can 
be computed as the two-dimensional nullspace of the linearisation of (A 2) in (x, 77 )oici. One 
component of this tangent space, spanned by ti = (x,?7)tan.b points along the circle. We 
denote its orthogonal complement in the tangent space by /2 = (x, 77 )tan,2, and take a small 
step (of size 5) from (x, 77)01^ along t2 to obtain an initial guess (x, 77) new. Newton iteration 
is used to solve (A 2), augmented with the two linear equations 

0=/ •^L,*:[-^<;'(0--^new(0]d^ + ^L,fcte~^°ew] for A; = 1, 2. (A4) 
Jo 

This gives for every (j) G [0, 2n] a new solution x^ on the resonance surface. The family 
of new solutions forms again a topological circle, so that the above procedure can be 
repeated to 'grow' the resonance surface as a family of topological circles. In practice one 
has to perform computations only for angles (j) G [0,27t/£], since for any x,^ the solution 
t X(i,{t — 77^)mod[o,i] ' j — 1 • • • ^ — 1, is also a solution. 
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